// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtsw4imgFileFormat.C                           //
// ************************************************************************* //

#include <avtsw4imgFileFormat.h>

#include <string>

#include <vtkDoubleArray.h>
#include <vtkFloatArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>

#include <DBOptionsAttributes.h>
#include <Expression.h>

#include <InvalidVariableException.h>
#include <InvalidDBTypeException.h>

#include <DebugStream.h>
#include <stdio.h>

#include <fcntl.h>
#include <visit-config.h>

using namespace std;

#ifdef _WIN32
#include <io.h>
#define OPEN  ::_open
#define CLOSE ::_close
#define LSEEK ::_lseek
#else
#include <unistd.h>
#define O_BINARY 0
#define OPEN  ::open
#define CLOSE ::close
#define LSEEK ::lseek
#endif

// ****************************************************************************
//  Method: avtsw4imgFileFormat constructor
//
//  Programmer: petersson1 -- generated by xml2avt
//  Creation:   Tue Nov 11 12:42:37 PDT 2014
//
// ****************************************************************************

avtsw4imgFileFormat::avtsw4imgFileFormat(const char *filename)
    : avtSTMDFileFormat(&filename, 1)
{
    // INITIALIZE DATA MEMBERS
    m_filename = filename;

    m_initialized = false;
    m_data_stored = false; // set to true if you want the solution field to be stored internally in this object (?)
}


// ****************************************************************************
//  Method: avtsw4imgFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: petersson1 -- generated by xml2avt
//  Creation:   Tue Nov 11 12:42:37 PDT 2014
//
// ****************************************************************************

void
avtsw4imgFileFormat::FreeUpResources(void)
{
}


// ****************************************************************************
//  Method: avtsw4imgFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: petersson1 -- generated by xml2avt
//  Creation:   Tue Nov 11 12:42:37 PDT 2014
//
// ****************************************************************************

void
avtsw4imgFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
{
    debug4 << "SW4> " << "Populatedatabase begins" << endl;

    if( !m_initialized )
        Initialize();

    debug4 << "SW4> " << "Populatedatabase after Initialize" << endl;
    // Mesh information
    avtMeshMetaData *mmd = new avtMeshMetaData;

    if (m_CartGrid)
    {
        mmd->name = "Cartesian grid";
    }
    else
    {
        mmd->name = "Curvilinear grid";
    }

    mmd->spatialDimension = 3;
    mmd->topologicalDimension = 3;

    mmd->numBlocks = m_nblocks;
        
    //   mmd->hasUnits = true;
    //   mmd->units = "meters";

    if (m_CartGrid)
    {
        mmd->meshType = AVT_RECTILINEAR_MESH;
    }
    else
    {
        mmd->meshType = AVT_CURVILINEAR_MESH;
    }

    debug4 << "SW4> " << "before b-loop " << endl;
    mmd->blockNames.resize(m_nblocks);
    for( int b=0 ; b < m_nblocks ; b++ )
    {
        debug4 << "SW4> " << "Block = " << b << endl;
        debug4 << "SW4> " << "grid size = " << m_gridsize[b] << endl;
        debug4 << "SW4> " << "no of pts = " << m_ni[b] << " " << m_nj[b]
               << " " << m_nk[b] << endl;
        char buf[50];
        snprintf(buf,50,"block=%i, h=%f", b, m_gridsize[b]);
        mmd->blockNames[b] = buf;
    }
    md->Add(mmd);
    debug4 << "SW4> " << "after b-loop " << endl;

    // Scalar variable information
    avtScalarMetaData *smd1 = new avtScalarMetaData;
    smd1->name = m_mode;
    // copy the mesh name from above
    smd1->meshName = mmd->name;
    smd1->centering = AVT_NODECENT;
    smd1->hasUnits = false;
    md->Add(smd1);
    debug4 << "SW4> " << "Populatedatabase done " << endl;
}


// ****************************************************************************
//  Method: avtsw4imgFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: petersson1 -- generated by xml2avt
//  Creation:   Tue Nov 11 12:42:37 PDT 2014
//
// ****************************************************************************

vtkDataSet *
avtsw4imgFileFormat::GetMesh(int domain, const char *meshname)
{
    if( !m_initialized )
        Initialize();

    debug4 << "SW4> " << "entering GetMesh " << endl;
    debug4 << "SW4> " << "meshname = '" << meshname << "'" << " domain = "
           << domain << endl;

    if (m_CartGrid || domain != m_nblocks-1)
    {
        debug4 << "SW4> " << "constructing a Cartesian grid, m_ni = "
               << m_ni[domain] << " m_nj = " << m_nj[domain] << " m_nk = "
               << m_nk[domain] << " m_gridsize = " << m_gridsize[domain]
               << endl;

        vtkDoubleArray *coords[3]= {0,0,0};
        coords[0] = vtkDoubleArray::New();
        coords[0]->SetNumberOfTuples(m_ni[domain]);
        coords[1] = vtkDoubleArray::New();
        coords[1]->SetNumberOfTuples(m_nj[domain]);
        coords[2] = vtkDoubleArray::New();
        coords[2]->SetNumberOfTuples(m_nk[domain]);
        // get the internal vtk data array pointers
        double* x = (double *)coords[0]->GetVoidPointer(0);
        double* y = (double *)coords[1]->GetVoidPointer(0);
        double* z = (double *)coords[2]->GetVoidPointer(0);
        // assign coordinates
        for( int i=0 ; i < m_ni[domain] ; i++ )
           x[i] = i*m_gridsize[domain]+m_xmin[domain];
        for( int i=0 ; i < m_nj[domain] ; i++ )
           y[i] = i*m_gridsize[domain]+m_ymin[domain];
        for( int i=0 ; i < m_nk[domain] ; i++ )
           z[i] = i*m_gridsize[domain]+m_zmin[domain];

        debug4 << "SW4> " << "done assigning x and y" << endl;
        
        vtkRectilinearGrid *rgrid= vtkRectilinearGrid::New();
        int dims[3];
        dims[0] = m_ni[domain];
        dims[1] = m_nj[domain];
        dims[2] = m_nk[domain];
        rgrid->SetDimensions(dims);
        rgrid->SetXCoordinates(coords[0]);
        rgrid->SetYCoordinates(coords[1]);
        rgrid->SetZCoordinates(coords[2]);
        coords[0]->Delete();
        coords[1]->Delete();
        coords[2]->Delete();
        debug4 << "SW4> " << "done making a Cartesian mesh " << endl;
        return rgrid;
    }
    else // curvilinear grid in Z for the last block
    {
        char errmsg[500];
        if( domain != m_nblocks-1 )
        {
            snprintf(errmsg, 500, "Error getMesh() domain = %i, m_nblocks = %i",
                domain, m_nblocks);
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }

        debug4 << "SW4> " << "constructing a curvilinear grid, m_CartGrid = "
               << m_CartGrid << " domain = " << domain << endl; 
        debug4 << "SW4> " << "m_ni = " << m_ni[domain] << " m_nj = "
               << m_nj[domain] << " m_nk = " << m_nk[domain] << endl;

        int dims[3] = {1,1,1};
        dims[0] = m_ni[domain];
        dims[1] = m_nj[domain];
        dims[2] = m_nk[domain];
        int nnodes = dims[0]*dims[1]*dims[2];
        debug4 << "SW4> " << "dims[0] = " << dims[0] << " dims[1] = " << dims[1] << " dims[2] = " << dims[2] << endl;

        // read z-coordinate from grid block on the sw4img file
        int fd = OPEN(m_filename.c_str(),O_RDONLY|O_BINARY);
        if( fd == -1 )
        {
            snprintf(errmsg,500,"Error opening grid file %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }

        // locate the grid block on file
        off_t nr = LSEEK(fd,m_gridOffset,SEEK_CUR);
        if( nr != m_gridOffset )
        {
            CLOSE(fd);
            snprintf(errmsg,500,"Error accessing grid array in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }

        // Create the structured grid.
        vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
        vtkPoints * points = vtkPoints::New();
        sgrid->SetPoints(points);
        sgrid->SetDimensions(dims);
        points->Delete();

        // get the data
        if( m_prec == 4 )
        {
            float *zarray = new float[nnodes];
            nr = read(fd,zarray,sizeof(float)*nnodes);
            if( (size_t)nr != sizeof(float)*nnodes )
            {
                CLOSE(fd);
                delete [] zarray;
                sgrid->Delete();
                snprintf(errmsg,500,"Error reading grid array in %s" , m_filename.c_str());
                EXCEPTION1( InvalidDBTypeException, errmsg );
            }

            points->SetDataType(VTK_FLOAT);
            points->SetNumberOfPoints(nnodes);

            float *pts = (float *) points->GetVoidPointer(0);
            float *zc = zarray;

            for(int k = 0; k < dims[2]; ++k)
            {
                for(int j = 0; j < dims[1]; ++j)
                {
                    for(int i = 0; i < dims[0]; ++i)
                    {
                        *pts++ = i*m_gridsize[domain]+m_xmin[domain]; // Cartesian in x and y
                        *pts++ = j*m_gridsize[domain]+m_ymin[domain];
                        *pts++ = *zc++;
                    }
                }
            }

            delete [] zarray;
        }
        else
        {
            double *zarray = new double[nnodes];
            nr = read(fd,zarray,sizeof(double)*nnodes);
            if( (size_t)nr != sizeof(double)*nnodes )
            {
                CLOSE(fd);
                delete [] zarray;
                sgrid->Delete();
                snprintf(errmsg,500,"Error reading dp grid array in %s" , m_filename.c_str());
                EXCEPTION1( InvalidDBTypeException, errmsg );
            }

            points->SetDataType(VTK_DOUBLE);
            points->SetNumberOfPoints(nnodes);

            double *pts = (double *) points->GetVoidPointer(0);
            double *zc = zarray;

            for(int k = 0; k < dims[2]; ++k)
            {
                for(int j = 0; j < dims[1]; ++j)
                {
                    for(int i = 0; i < dims[0]; ++i)
                    {
                        *pts++ = i*m_gridsize[domain]+m_xmin[domain]; // Cartesian in x and y
                        *pts++ = j*m_gridsize[domain]+m_ymin[domain];
                        *pts++ = *zc++;
                    }
                }
            }

            delete [] zarray;
        }

        // end reading z-ccordinates: close the file
        CLOSE(fd);

        debug4 << "SW4> " << "done making a curvilinear mesh" << endl;
        return sgrid;
    }
}


// ****************************************************************************
//  Method: avtsw4imgFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: petersson1 -- generated by xml2avt
//  Creation:   Tue Nov 11 12:42:37 PDT 2014
//
// ****************************************************************************

vtkDataArray *
avtsw4imgFileFormat::GetVar(int domain, const char *varname)
{
    debug4 << "SW4> " << "entering GetVar " << endl;
    if( !m_initialized )
        Initialize();

    debug4 << "SW4> " << "after Initialize() " << endl;
    char errmsg[500];

    vtkDataArray* ret = NULL;

    int fd = OPEN(m_filename.c_str(),O_RDONLY|O_BINARY);
    if( fd == -1 )
    {
        snprintf(errmsg,500,"Error opening file %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    off_t nr = LSEEK(fd,m_offset[domain],SEEK_CUR);
    if( nr != m_offset[domain] )
    {
        CLOSE(fd);
        snprintf(errmsg,500,"Error accessing array in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    if( m_prec == 4 )
    {
        vtkFloatArray* arr = vtkFloatArray::New();
        size_t npts = (size_t) m_ni[domain]*m_nj[domain]*m_nk[domain];
        arr->SetNumberOfTuples(npts);
        float* data = (float *)arr->GetVoidPointer(0);

        nr = read(fd,data,sizeof(float)*npts);
        if( (size_t)nr != sizeof(float)*npts )
        {
            CLOSE(fd);
            arr->Delete();
            snprintf(errmsg,500,"Error reading array in %s" , m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }

        ret = arr;
    }
    else
    {
        vtkDoubleArray* arr = vtkDoubleArray::New();
        size_t npts = (size_t) m_ni[domain]*m_nj[domain]*m_nk[domain];
        arr->SetNumberOfTuples(npts);
        double* data = (double *)arr->GetVoidPointer(0);

        nr = read(fd,data,sizeof(double)*npts);
        if( (size_t)nr != sizeof(double)*npts )
        {
            CLOSE(fd); 
            arr->Delete();
            snprintf(errmsg,500,"Error reading array in %s" , m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }

        ret = arr;
    }

    // close the file
    CLOSE(fd);

    debug4 << "SW4> " << "leaving GetVar " << endl;
    return ret;
}


// ****************************************************************************
//  Method: avtsw4imgFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: petersson1 -- generated by xml2avt
//  Creation:   Tue Nov 11 12:42:37 PDT 2014
//
// ****************************************************************************

vtkDataArray *
avtsw4imgFileFormat::GetVectorVar(int domain, const char *varname)
{
    // there are no vector variables
    return 0;
}

void avtsw4imgFileFormat::Initialize()
{
    debug4 << "SW4> " << "entering Initialize()" << endl;
    char errmsg[500];

    // should check for 3D in the filename and handle that case within this reader
    size_t k=m_filename.rfind('.');
    size_t find3D = m_filename.find("3Dimg");
    if (find3D != string::npos && find3D == k+1)
    {
        // We have a 'mode.3Dimg' type file -- this is for volimage
        debug4 << "SW4> " << "opening volimage file, continuing..." << endl;
        m_volimage = true;
    }
    else
    {
        debug4 << "SW4> " << "Not a volimage file, continuing..." << endl; 
        m_volimage = false;
    }

    int fd = OPEN( m_filename.c_str(), O_RDONLY|O_BINARY );
    if( fd == -1 )
    {
        snprintf(errmsg,500,"Error opening file %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    debug4 << "SW4> " << "file opened " << endl;

    // reading the header...
    // precision
    size_t nr = read(fd,&m_prec,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading precision in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    // number of patches (blocks)
    nr = read(fd,&m_nblocks,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading nblocks in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    // time
    nr = read(fd,&m_time,sizeof(double) );
    if( nr != sizeof(double) )
    {
        snprintf(errmsg,500,"Error reading time in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    debug4 << "SW4> " << "prec = " << m_prec << " " << " m_nblocks = " << m_nblocks << " time= " << m_time << endl;
    // plane (0=x-const, 1=y-const, 2=z-const)
    nr = read(fd,&m_plane,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading plane in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    // plane_coord
    nr = read(fd,&m_plane_coord,sizeof(double) );
    if( nr != sizeof(double) )
    {
        snprintf(errmsg,500,"Error reading plane_coordinate in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    // mode
    nr = read(fd,&m_mode_number,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading mode number in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    switch (m_mode_number)
    {
        case 1:
            m_mode="ux"; break;
        case 2:
            m_mode="uy"; break;
        case 3:
            m_mode="uz"; break;
        case 4:
            m_mode="rho"; break;
        case 5:
            m_mode="lambda"; break;
        case 6:
            m_mode="mu"; break;
        case 7:
            m_mode="cp"; break;
        case 8:
            m_mode="cs"; break;
        case 9:
            m_mode="uex"; break;
        case 10:
            m_mode="uey"; break;
        case 11:
            m_mode="uez"; break;
        case 12:
            m_mode="div"; break;
        case 13:
            m_mode="curl"; break;
        case 14:
            m_mode="divt"; break;
        case 15:
            m_mode="curlt"; break;
        case 16:
            m_mode="lat"; break;
        case 17:
            m_mode="lon"; break;
        case 18:
            m_mode="topo"; break;
        case 19:
            m_mode="x"; break;
        case 20:
            m_mode="y"; break;
        case 21:
            m_mode="z"; break;
        case 22:
            m_mode="uxerr"; break;
        case 23:
            m_mode="uyerr"; break;
        case 24:
            m_mode="uzerr"; break;
        case 25:
            m_mode="magt"; break;
        case 26:
            m_mode="hmagt"; break;
        case 27:
            m_mode="hmaxt"; break;
        case 28:
            m_mode="vmaxt"; break;
        case 29:
            m_mode="mag"; break;
        case 30:
            m_mode="hmag"; break;
        case 31:
            m_mode="hmax"; break;
        case 32:
            m_mode="vmax"; break;

        default:
            m_mode="unknown"; break;
    }
    debug4 << "SW4> " << "mode number = " << m_mode_number << " mode = "
           << m_mode << endl;
    
    // grdinfo
    nr = read(fd,&m_grdinfo,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading grdinfo in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    m_CartGrid = (m_grdinfo == 0);
    
    debug4 << "SW4> " << "grid info = " << m_grdinfo << endl;

    // file creation time    
    char cgfn[25];
    nr = read(fd,&cgfn,25*sizeof(char) );
    if( nr != 25*sizeof(char) )
    {
       snprintf(errmsg, 500, "Error reading creation time stamp in %s",
           m_filename.c_str());
       EXCEPTION1( InvalidDBTypeException, errmsg );
    }

    // the c-string is terminated by a null character
    cgfn[24]='\0';
    m_file_creation_time = cgfn;
    debug4 << "SW4> " << "reading creation time stamp='"
           << m_file_creation_time << "'" << endl;

    // end of first part of header

    // 2nd part of header: grid size
    m_gridsize.resize(m_nblocks);
    m_xmin.resize(m_nblocks);
    m_ymin.resize(m_nblocks);
    m_zmin.resize(m_nblocks);
    m_z0.resize(m_nblocks);
    m_ni.resize(m_nblocks);
    m_nj.resize(m_nblocks);
    m_nk.resize(m_nblocks);
    m_offset.resize(m_nblocks);

    // size of the first part of the header
    off_t header_offset = 5*sizeof(int)+2*sizeof(double)+25*sizeof(char);
    int dims[6];
    debug4 << "SW4> " << "reading block headers " << endl;
    for( int b=0 ; b < m_nblocks ; b++ )
    {
        debug4 << "SW4> " << "b = " << b << endl;
        // grid size
        nr = read(fd,&m_gridsize[b],sizeof(double));
        if( nr != sizeof(double) )
        {
            snprintf(errmsg,500,"Error reading gridsizes in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        // z0
        nr = read(fd,&m_z0[b],sizeof(double));
        if( nr != sizeof(double) )
        {
            snprintf(errmsg,500,"Error reading z0 in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        int expect_dim = (m_volimage) ? 6 : 4;
        nr = read(fd,dims,sizeof(int)*expect_dim);
        if( nr != sizeof(int)*expect_dim )
        {
            snprintf(errmsg,500,"Error reading dimensions in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        if( m_plane == 0 )
        {
            m_ni[b] = 1;
            m_nj[b] = dims[1]-dims[0]+1;
            m_nk[b] = dims[3]-dims[2]+1;
            m_xmin[b] = m_plane_coord;
            m_ymin[b] = (dims[0]-1)*m_gridsize[b];
            m_zmin[b] = m_z0[b] + (dims[2]-1)*m_gridsize[b];
        }
        else if( m_plane == 1 )
        {
            m_ni[b] = dims[1]-dims[0]+1;
            m_nj[b] = 1;
            m_nk[b] = dims[3]-dims[2]+1;
            m_xmin[b] = (dims[0]-1)*m_gridsize[b];
            m_ymin[b] = m_plane_coord;
            m_zmin[b] = m_z0[b] + (dims[2]-1)*m_gridsize[b];
        }
        else
        {
            m_ni[b] = dims[1]-dims[0]+1;
            m_nj[b] = dims[3]-dims[2]+1;
            m_nk[b] = (m_volimage) ? dims[5]-dims[4]+1 : 1;
            m_xmin[b] = (dims[0]-1)*m_gridsize[b];
            m_ymin[b] = (dims[2]-1)*m_gridsize[b];
            m_zmin[b] = (m_volimage) ?
                m_z0[b] + (dims[4]-1)*m_gridsize[b] : m_plane_coord;
        }
        debug4 << "SW4> " << "b = " << b << " dims[0] = " << dims[0]
               << " dims[1] = " << dims[1] << " dims[2] = " << dims[2]
               << " dims[3] = " << dims[3];
        if (m_volimage)
        {
            debug4 << " dims[4] = " << dims[4] << " dims[5] = " << dims[5];
        }
        debug4 << endl;

        debug4 << "SW4> " << "z0 = " << m_z0[b] << " x_min = "
               << m_xmin[b] << " y_min = " << m_ymin[b]
               << " z_min = " << m_zmin[b] << endl;
        // old
        header_offset += 2*sizeof(double) + expect_dim*sizeof(int);
        debug4 << "SW4> " << "block = " << b << " dims[0] = " << dims[0]
               << " dims[1] = " << dims[1] << " dims[2] = "
               << dims[2] << " dims[3] = " << dims[3];
        if (m_volimage)
        {
            debug4 << " dims[4] = " << dims[4] << " dims[5] = " << dims[5];
        }
        debug4 << endl;
        debug4 << "SW4> " << " z_min = " << m_zmin[b] << endl;
    }
    // end 2nd part of header

    m_offset[0] = header_offset;
    int datasize;
    if( m_prec == 4 )
        datasize = sizeof(float);
    else
        datasize = sizeof(double);

    // calculate offsets for accessing data
    for( int b=1 ; b < m_nblocks ; b++ )
        m_offset[b] = m_offset[b-1] + ((size_t) m_ni[b-1])*m_nj[b-1]*m_nk[b-1]*datasize;

    // assign m_gridOffset
    if (m_CartGrid)
    {
        m_gridOffset = -1;
    }
    else
    {
        // lb : index of the last block
        int lb = m_nblocks-1;
        m_gridOffset = m_offset[lb] + ((size_t) m_ni[lb])*m_nj[lb]*m_nk[lb]*datasize;
    }

    debug4 << "SW4> " << " done headers, data stored =  " << m_data_stored << endl;

    CLOSE(fd); // closing the solution file

    // done intializing...
    m_initialized = true;
    debug4 << "SW4> " << "leaving Initialize()" << endl;
}
